This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

1. Variable parameters in theta_matrix

1.1

#
environ_params <- list(
  M = 6,        ## Actors
  N = 10,       ## Components
  BI_PROB = 0, ## Environmental Density (DGP hyperparameter)
  component_matrix_start = 'rand', ##**TODO** Implement: 'rand','modular','semi-modular',...
  rand_seed = 1234,
  plot_init = F,
  name = '_test_tutorial_nb_'
)
#
env3 <- SaomNkRSienaBiEnv$new(environ_params)
## 
## TEST FROM CALLED CLASS INIT *BEFORE* BASE INIT
## 
## CALLED _BASE_ INIT
## 
## TEST FROM CALLED CLASS INIT *AFTER* BASE INIT
#max_waves <- 15
#array_bi_net <-  env3$bi_env_arr[ , , 1:min(max_waves, dim(env3$bi_env_arr)[3]) ] ##**MAIN STEP**
# self$bi_env_arr[,,1] + self$bi_env_arr[,,2] + self$bi_env_arr[,,3] 
# self$bi_env_arr[,,4] + self$bi_env_arr[,,5] + self$bi_env_arr[,,6] 
array_bi_net <- array(c(env3$bipartite_matrix, env3$bipartite_matrix), 
                      dim=c(env3$M,  env3$N, 2))
print(array_bi_net)
## , , 1
## 
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    0    0    0    0    0    0    0    0    0     0
## [2,]    0    0    0    0    0    0    0    0    0     0
## [3,]    0    0    0    0    0    0    0    0    0     0
## [4,]    0    0    0    0    0    0    0    0    0     0
## [5,]    0    0    0    0    0    0    0    0    0     0
## [6,]    0    0    0    0    0    0    0    0    0     0
## 
## , , 2
## 
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    0    0    0    0    0    0    0    0    0     0
## [2,]    0    0    0    0    0    0    0    0    0     0
## [3,]    0    0    0    0    0    0    0    0    0     0
## [4,]    0    0    0    0    0    0    0    0    0     0
## [5,]    0    0    0    0    0    0    0    0    0     0
## [6,]    0    0    0    0    0    0    0    0    0     0
strategies <- list(
  egoX   =   c(-1, 0, 1),
  inPopX =   c( 1, 0, -1)
)


## 2.b. Component Payoffs vector
set.seed(12345)
component_payoffs <-  runif(environ_params$N, min = 0, max = 1)
## 2. Strategies sets the objective function as a linear combination of network stats across DVs
#
actor_strats_list <- lapply(strategies, function(strat) rep(strat,  environ_params$M/length(strat)) )

#
structure_model <- list(
  dv_bipartite = list(
    name = 'self$bipartite_rsienaDV',
    effects = list( ##**STRUCTURAL EFFECTS -- dyadic/network endogeneity sources**
      list(effect='density', parameter= 0, fix=T, dv_name=DV_NAME), ##interaction1 = NULL
      list(effect='inPop',   parameter= 0, fix=T, dv_name=DV_NAME), #interaction1 = NUL
      list(effect='outAct',  parameter= 0, fix=T, dv_name=DV_NAME)
    ),
    ## COVARIATE EFFECTS
    coCovars = list( 
      ##** COMPONENTS : MONADIC CONSTANT COVARIATE EFFECTS **##
      list(effect='altX',   parameter= .1,  fix=T,dv_name=DV_NAME, 
           interaction1='self$component_1_coCovar', x = component_payoffs 
      ),
      ##** STRATEGIES : MONADIC CONSTANT COVARIATE EFFECTS **##
      list(effect='egoX',   parameter= .01,  fix=T,dv_name=DV_NAME, 
           interaction1='self$strat_1_coCovar',     x = actor_strats_list[[1]] 
      ), #interaction1 = NULL
      list(effect='inPopX', parameter= .1,  fix=T,dv_name=DV_NAME, 
           interaction1='self$strat_2_coCovar',  x = actor_strats_list[[2]] 
      )
    ),
    varCovars = list() ##**MONADIC TIME-VARYING COVARIATE EFFECTS -- DYNAMIC STRATEGY PROGRAMS**
  )
)

#
structeffs <- structure_model$dv_bipartite$effects %>% ldply(as.data.frame) # use previous effects
coCovars   <- structure_model$dv_bipartite$coCovars %>% ldply(function(x)as.data.frame(x[! names(x)%in%c('interaction1','x')])) 
effects    <- structeffs %>% bind_rows(coCovars)
# names_theta_in <- effects$effectName[effects$include][params_norates]
theta_in       <- effects$parameter

## Number of decision chain steps to simulate
nrows <- environ_params$M * 10
#
ncols <- length(theta_in)
shock.chunks <- 4  ## before|after
theta_matrix <- matrix(NA,  nrow = nrows, ncol=ncols  )
for (j in 1:ncols) {
  nchunks   <- floor( nrows / shock.chunks ) 
  beta_seq_1 <- rep(  1 * theta_in[j], nchunks * 1 )  ## How many pd_chunks at what theta value?
  beta_seq_2 <- rep(  1 * theta_in[j], nchunks * 1 ) ## number of pd_chunks at what theta value?
  beta_seq_3 <- rep(  1 * theta_in[j], nchunks * 2 ) ## number of pd_chunks at what theta value?
  vec <- c( beta_seq_1, beta_seq_2, beta_seq_3 )
  rowdiff <- nrow(theta_matrix) - length(vec)
  if (rowdiff != 0) {
    vec <- c( vec, rep(vec[length(vec)], round(abs(rowdiff))) ) ## add the last elements again if too short
  }
  #
  theta_matrix[, j] <- vec
}

print(theta_matrix)
##       [,1] [,2] [,3] [,4] [,5] [,6]
##  [1,]    0    0    0  0.1 0.01  0.1
##  [2,]    0    0    0  0.1 0.01  0.1
##  [3,]    0    0    0  0.1 0.01  0.1
##  [4,]    0    0    0  0.1 0.01  0.1
##  [5,]    0    0    0  0.1 0.01  0.1
##  [6,]    0    0    0  0.1 0.01  0.1
##  [7,]    0    0    0  0.1 0.01  0.1
##  [8,]    0    0    0  0.1 0.01  0.1
##  [9,]    0    0    0  0.1 0.01  0.1
## [10,]    0    0    0  0.1 0.01  0.1
## [11,]    0    0    0  0.1 0.01  0.1
## [12,]    0    0    0  0.1 0.01  0.1
## [13,]    0    0    0  0.1 0.01  0.1
## [14,]    0    0    0  0.1 0.01  0.1
## [15,]    0    0    0  0.1 0.01  0.1
## [16,]    0    0    0  0.1 0.01  0.1
## [17,]    0    0    0  0.1 0.01  0.1
## [18,]    0    0    0  0.1 0.01  0.1
## [19,]    0    0    0  0.1 0.01  0.1
## [20,]    0    0    0  0.1 0.01  0.1
## [21,]    0    0    0  0.1 0.01  0.1
## [22,]    0    0    0  0.1 0.01  0.1
## [23,]    0    0    0  0.1 0.01  0.1
## [24,]    0    0    0  0.1 0.01  0.1
## [25,]    0    0    0  0.1 0.01  0.1
## [26,]    0    0    0  0.1 0.01  0.1
## [27,]    0    0    0  0.1 0.01  0.1
## [28,]    0    0    0  0.1 0.01  0.1
## [29,]    0    0    0  0.1 0.01  0.1
## [30,]    0    0    0  0.1 0.01  0.1
## [31,]    0    0    0  0.1 0.01  0.1
## [32,]    0    0    0  0.1 0.01  0.1
## [33,]    0    0    0  0.1 0.01  0.1
## [34,]    0    0    0  0.1 0.01  0.1
## [35,]    0    0    0  0.1 0.01  0.1
## [36,]    0    0    0  0.1 0.01  0.1
## [37,]    0    0    0  0.1 0.01  0.1
## [38,]    0    0    0  0.1 0.01  0.1
## [39,]    0    0    0  0.1 0.01  0.1
## [40,]    0    0    0  0.1 0.01  0.1
## [41,]    0    0    0  0.1 0.01  0.1
## [42,]    0    0    0  0.1 0.01  0.1
## [43,]    0    0    0  0.1 0.01  0.1
## [44,]    0    0    0  0.1 0.01  0.1
## [45,]    0    0    0  0.1 0.01  0.1
## [46,]    0    0    0  0.1 0.01  0.1
## [47,]    0    0    0  0.1 0.01  0.1
## [48,]    0    0    0  0.1 0.01  0.1
## [49,]    0    0    0  0.1 0.01  0.1
## [50,]    0    0    0  0.1 0.01  0.1
## [51,]    0    0    0  0.1 0.01  0.1
## [52,]    0    0    0  0.1 0.01  0.1
## [53,]    0    0    0  0.1 0.01  0.1
## [54,]    0    0    0  0.1 0.01  0.1
## [55,]    0    0    0  0.1 0.01  0.1
## [56,]    0    0    0  0.1 0.01  0.1
## [57,]    0    0    0  0.1 0.01  0.1
## [58,]    0    0    0  0.1 0.01  0.1
## [59,]    0    0    0  0.1 0.01  0.1
## [60,]    0    0    0  0.1 0.01  0.1
## Run Rsiena search using variable parameters in theta_matrix
env3$search_rsiena_shocks(
  array_bi_net,
  structure_model,
  theta_matrix,
  digits = 4,
  run_seed = 12345
)
## [1] "self$rsiena_data : "
## Dependent variables:  self$bipartite_rsienaDV 
## Number of observations: 2 
## 
## Nodesets                 ACTORS      COMPONENTS 
## Number of nodes               6              10 
## 
## Dependent variable self$bipartite_rsienaDV
## Type               bipartite              
## Observations       2                      
## First nodeset      ACTORS                 
## Second nodeset     COMPONENTS             
## Densities          NA NA                  
## 
## Constant covariates:  self$component_1_coCovar, self$strat_1_coCovar, self$strat_2_coCovar 
## 
##  structural effects i=1, j=1
## $effect
## [1] "density"
## 
## $parameter
## [1] 0
## 
## $fix
## [1] TRUE
## 
## $dv_name
## [1] "self$bipartite_rsienaDV"
## 
##   effectName          include fix  test  initialValue parm
## 1 outdegree (density) TRUE    TRUE FALSE   -1.60944   0   
##   effectName          include fix  test  initialValue parm
## 1 outdegree (density) TRUE    TRUE FALSE          0   0   
## 
##  structural effects i=1, j=2
## $effect
## [1] "inPop"
## 
## $parameter
## [1] 0
## 
## $fix
## [1] TRUE
## 
## $dv_name
## [1] "self$bipartite_rsienaDV"
## 
##   effectName            include fix  test  initialValue parm
## 1 indegree - popularity TRUE    TRUE FALSE          0   0   
##   effectName            include fix  test  initialValue parm
## 1 indegree - popularity TRUE    TRUE FALSE          0   0   
## 
##  structural effects i=1, j=3
## $effect
## [1] "outAct"
## 
## $parameter
## [1] 0
## 
## $fix
## [1] TRUE
## 
## $dv_name
## [1] "self$bipartite_rsienaDV"
## 
##   effectName           include fix  test  initialValue parm
## 1 outdegree - activity TRUE    TRUE FALSE          0   0   
##   effectName           include fix  test  initialValue parm
## 1 outdegree - activity TRUE    TRUE FALSE          0   0   
## 
##  coCovars i=1, j=1
## $effect
## [1] "altX"
## 
## $parameter
## [1] 0.1
## 
## $fix
## [1] TRUE
## 
## $dv_name
## [1] "self$bipartite_rsienaDV"
## 
## $interaction1
## [1] "self$component_1_coCovar"
## 
## $x
##  [1] 0.7209039 0.8757732 0.7609823 0.8861246 0.4564810 0.1663718 0.3250954
##  [8] 0.5092243 0.7277053 0.9897369
## 
##   effectName                     include fix  test  initialValue parm
## 1 self$component_1_coCovar alter TRUE    TRUE FALSE          0   0   
##   effectName                     include fix  test  initialValue parm
## 1 self$component_1_coCovar alter TRUE    TRUE FALSE          0   0.1 
## 
##  coCovars i=1, j=2
## $effect
## [1] "egoX"
## 
## $parameter
## [1] 0.01
## 
## $fix
## [1] TRUE
## 
## $dv_name
## [1] "self$bipartite_rsienaDV"
## 
## $interaction1
## [1] "self$strat_1_coCovar"
## 
## $x
## [1] -1  0  1 -1  0  1
## 
##   effectName               include fix  test  initialValue parm
## 1 self$strat_1_coCovar ego TRUE    TRUE FALSE          0   0   
##   effectName               include fix  test  initialValue parm
## 1 self$strat_1_coCovar ego TRUE    TRUE FALSE          0   0.01
## 
##  coCovars i=1, j=3
## $effect
## [1] "inPopX"
## 
## $parameter
## [1] 0.1
## 
## $fix
## [1] TRUE
## 
## $dv_name
## [1] "self$bipartite_rsienaDV"
## 
## $interaction1
## [1] "self$strat_2_coCovar"
## 
## $x
## [1]  1  0 -1  1  0 -1
## 
##   effectName                                    include fix  test  initialValue
## 1 ind. pop.^(1/#) weighted self$strat_2_coCovar TRUE    TRUE FALSE          0  
##   parm
## 1 1   
##   effectName                                    include fix  test  initialValue
## 1 ind. pop.^(1/#) weighted self$strat_2_coCovar TRUE    TRUE FALSE          0  
##   parm
## 1 0.1 
## If you use this algorithm object, siena07 will create/use an output file _test_tutorial_nb__173793437558.txt .
## 
## Start phase 0 
## theta: 0 0 0 0 0 0 
## 
## Start phase 3 
## Parameter values used for simulations
## 
##                                                           Mean      Standard      
##                                                             value   Deviation     
## 
## Rate parameters: 
##   0       Rate parameter                                    NA    ( NA        )   
## 
## Other parameters: 
##   1. eval outdegree (density)                             0.00    (  0        )   
##   2. eval indegree - popularity                           0.00    (  0        )   
##   3. eval outdegree - activity                            0.00    (  0        )   
##   4. eval self$component_1_coCovar alter                  0.10    (  0        )   
##   5. eval self$strat_1_coCovar ego                        0.01    (  0        )   
##   6. eval ind. pop.^(1/0.1) weighted self$strat_2_coCovar 0.10    (  0        )   
## 
## Simulated means and standard deviations
##   1. Number of ties                                             0.817    0.390 
##   2. Sum of squared indegrees                                   0.817    0.390 
##   3. Sum of squared outdegrees                                  0.817    0.390 
##   4. Sum of indegrees x self$component_1_coCovar               -0.002    0.199 
##   5. Sum of outdegrees x self$strat_1_coCovar                   0.017    0.725 
##   6. indegree pop.^(1/0.1) weighted self$strat_2_coCovar        0.000    0.611 
## 
## 
## Simulated statistics are in x$sf
## and simulated dependent variables in x$sims, where x is the created object.
## 
## Total of 60 iteration steps.
## 
## Covariance matrix of estimates (correlations below diagonal)
## 
##            0            0            0            0            0            0
##          NaN            0            0            0            0            0
##          NaN          NaN            0            0            0            0
##          NaN          NaN          NaN            0            0            0
##          NaN          NaN          NaN          NaN            0            0
##          NaN          NaN          NaN          NaN          NaN            0
## 
## Derivative matrix of expected statistics X by parameters:
## 
##        0.150        0.150        0.150        0.000        0.003        0.000
##        0.150        0.150        0.150        0.000        0.003        0.000
##        0.150        0.150        0.150        0.000        0.003        0.000
##        0.000        0.000        0.000        0.039        0.014       -0.013
##        0.000        0.000        0.000        0.001        0.047       -0.033
##        0.000        0.000        0.000        0.000        0.000        0.000
## 
## Covariance matrix of X (correlations below diagonal):
## 
##        0.152        0.152        0.152        0.000        0.003        0.000
##        1.000        0.152        0.152        0.000        0.003        0.000
##        1.000        1.000        0.152        0.000        0.003        0.000
##       -0.005       -0.005       -0.005        0.040        0.015       -0.013
##        0.011        0.011        0.011        0.102        0.525       -0.373
##        0.000        0.000        0.000       -0.109       -0.843        0.373
## 
## =====================================================================
##                                                  Model 1             
## ---------------------------------------------------------------------
## Rate parameter period 1                           0.1361 (0.1181)    
## outdegree (density)                               0.0000 (0.0000)    
## indegree - popularity                             0.0000 (0.0000)    
## outdegree - activity                              0.0000 (0.0000)    
## self$component_1_coCovar alter                    0.1000 (0.0000) ***
## self$strat_1_coCovar ego                          0.0100 (0.0000) ***
## ind. pop.^(1/0.1) weighted self$strat_2_coCovar   0.1000 (0.0000) ***
## ---------------------------------------------------------------------
## Iterations                                       60                  
## =====================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05
##
env3$search_rsiena_process_ministep_chain()
##   dv_type dv_type_bin              dv_varname id_from id_to beh_difference
## 1 Network           0 self$bipartite_rsienaDV       6     9              0
## 2 Network           0 self$bipartite_rsienaDV       2     6              0
## 3 Network           0 self$bipartite_rsienaDV       6     1              0
## 4 Network           0 self$bipartite_rsienaDV       5     1              0
## 5 Network           0 self$bipartite_rsienaDV       3     5              0
## 6 Network           0 self$bipartite_rsienaDV       2    11              0
##   reciprocal_rate LogOptionSetProb LogChoiceProb diagonal stability tie_change
## 1       0.1666667        -1.791759     -2.388695    FALSE     FALSE       TRUE
## 2       0.1666667        -1.791759     -2.445733    FALSE     FALSE       TRUE
## 3       0.1666667        -1.791759     -2.389375    FALSE     FALSE       TRUE
## 4       0.1666667        -1.791759     -2.390279    FALSE     FALSE       TRUE
## 5       0.1666667        -1.791759     -2.415817    FALSE     FALSE       TRUE
## 6       0.1666667        -1.791759     -2.398186    FALSE      TRUE      FALSE
##   chain_step_id
## 1             1
## 2             2
## 3             3
## 4             4
## 5             5
## 6             5
##    dv_type           dv_type_bin  dv_varname           id_from     
##  Length:60          Min.   :0    Length:60          Min.   :1.000  
##  Class :character   1st Qu.:0    Class :character   1st Qu.:2.000  
##  Mode  :character   Median :0    Mode  :character   Median :4.000  
##                     Mean   :0                       Mean   :3.667  
##                     3rd Qu.:0                       3rd Qu.:5.000  
##                     Max.   :0                       Max.   :6.000  
##      id_to       beh_difference reciprocal_rate  LogOptionSetProb
##  Min.   : 1.00   Min.   :0      Min.   :0.1667   Min.   :-1.792  
##  1st Qu.: 3.75   1st Qu.:0      1st Qu.:0.1667   1st Qu.:-1.792  
##  Median : 6.50   Median :0      Median :0.1667   Median :-1.792  
##  Mean   : 6.50   Mean   :0      Mean   :0.1667   Mean   :-1.792  
##  3rd Qu.: 9.00   3rd Qu.:0      3rd Qu.:0.1667   3rd Qu.:-1.792  
##  Max.   :11.00   Max.   :0      Max.   :0.1667   Max.   :-1.792  
##  LogChoiceProb      diagonal         stability       tie_change     
##  Min.   :-2.447   Length:60          Mode :logical   Mode :logical  
##  1st Qu.:-2.412   Class :character   FALSE:49        FALSE:11       
##  Median :-2.390   Mode  :character   TRUE :11        TRUE :49       
##  Mean   :-2.398                                                     
##  3rd Qu.:-2.387                                                     
##  Max.   :-2.362                                                     
##  chain_step_id  
##  Min.   : 1.00  
##  1st Qu.:15.75  
##  Median :30.50  
##  Mean   :30.30  
##  3rd Qu.:45.00  
##  Max.   :60.00  
## [1] 60 13
##
env3$search_rsiena_process_stats()

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

snapshot_ids <- round( seq(1, dim(env3$bi_env_arr)[3], length=3 ) )
for (i in 1:length(snapshot_ids)) {
  step <- snapshot_ids[ i ]
  env3$plot_bipartite_system_from_mat(env3$bi_env_arr[,,step], step)
}

# plotlist <- list()
# for (i in 1:length(snapshot_ids)) {
#   chain_step <- snapshot_ids[ i ]
#   plotlist[[i]] <-  self$plot_bipartite_system_from_mat(
#     self$bi_env_arr[,,chain_step], 
#     chain_step
#   ) 
# }
# ggarrange(plotlist = plotlist, nrow=length(plotlist))
## [1] 360   5
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

## [1] 360   5
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

##**TODO value_weighted** 
##**Multiply by the theta value in the function processing actor_stats_df **
#
env3$plot_utility_components(loess_span=0.35)
## `geom_smooth()` using formula = 'y ~ x'

## ggtitle('Actor Utility Contributed by Components')
##
env3$plot_utility_components(use_contributions=TRUE, loess_span=0.35)
## `geom_smooth()` using formula = 'y ~ x'

##
env3$plot_actor_degrees(loess_span = 0.25)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

##
env3$plot_component_degrees(loess_span = 0.35)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

######################
##**TODO**

self <- env3


#
config_param_vals <- c(
  unlist(sapply(self$config_structure_model$dv_bipartite$effects, function(x) x$parameter)),
  unlist(sapply(self$config_structure_model$dv_bipartite$coCovars, function(x) x$parameter)),
  unlist(sapply(self$config_structure_model$dv_bipartite$varCovars, function(x) x$parameter))
)
fixed_params <- c(
  unlist(sapply(self$config_structure_model$dv_bipartite$effects, function(x) x$fix)),
  unlist(sapply(self$config_structure_model$dv_bipartite$coCovars, function(x) x$fix)),
  unlist(sapply(self$config_structure_model$dv_bipartite$varCovars, function(x) x$fix))
)
#
#UTILITY
theta   <- self$rsiena_model$thetaUsed
## replace fixed theta param with given param values (instead of zero default when effect is fixed)
# if (any(fixed_params)) theta[ fixed_params ] <- config_param_vals[ fixed_params ]



## PLOT ACTOR UTILITY LANDSCAPE
utilist <- list()
statsl <- list()
pltlist <- list()

# step_ids <- round(c(1,8,64,512, self$rsiena_model$n3/2, self$rsiena_model$n3))
acpds <- seq(1, nrow(theta), by= self$M )
step_ids <- c(acpds[1:6])
for (step_id in step_ids) {
  # step_id <- 1
  cat(sprintf('\nstep %s: actors ', step_id))
  
  bi_env_mat_step <- self$bi_env_arr[, , step_id]
  
  act_counterfacts <- list()
  

  for (i in 1:self$M) {
    cat(sprintf(' %s ', i))
    
    
    ##**ACTOR i DECISION PERSPECTIVE** 0000000000000000000000000000000000
    
    ##**TODO**
    ## ALL actor-component counterfactual configuations for Actor i  (2^N rows by N cols)
    iland <- expand.grid(lapply(1:self$N, function(x) 0:1 ))
    iland_config_step_row_id <- which(apply(iland, 1, function(x) all(x == bi_env_mat_step[i,]) ))
    
    tmpmat <- bi_env_mat_step
    ## ifit dimensions [ M, 2^N ]
    ifit <- apply(iland, 1, function(x){
      tmpmat[i,] <- x  ## set counterfactual actor-component configuration
      config_fit <- if('matrix' %in% class(theta)) {
        self$get_struct_mod_stats_mat_from_bi_mat( tmpmat ) %*% theta[step_id,]
      } else {
        self$get_struct_mod_stats_mat_from_bi_mat( tmpmat ) %*% theta
      }
      return( config_fit ) ## compute utility vector
    }) 
    
    # apply
    
    ids.max <- which(ifit == max(ifit), arr.ind = TRUE) 
    fits.max <- ifit[ ids.max ]
    nmax <- length(fits.max)
    
    # ifit
    
    # ids.max[1,]
    
    act_counterfacts[[i]] <- list(
      ifit = ifit, #Given all other ties, Actor i's configurations applied to utility func for all actors
      ids.max = ids.max,
      fits.max = fits.max,
      nmax = nmax
    )
    
    ## distances of each counterfactual configuration 
    dist_counterfac <- iland - bi_env_mat_step[i, ]
    z <- apply( dist_counterfac, 1, function(x) {
      c(dist=sum(x!=0),`drop`=sum(x==-1), `nochange`=sum(x==0), `add`=sum(x==1), 
        counterfac=paste(x, collapse = '|'), start=paste(bi_env_mat_step[i, ], collapse = '|') )
    })
    # z <- z[,order(z['d',], decreasing = F)]
    ##**TODO**
    ##**INFORMATION LOSS**
    ##  This uses only ego's own counterfactual fits (landscape)
    ##  but the corresponding other actor's affected fits (landscapes are not currently used )
    z <- rbind(z, utility_ego=ifit[ i , ] )
    z <- rbind(z, utility_alter_mean= colMeans(ifit[ -i , ], na.rm = T) )
    z <- rbind(z, utility_alter_sd  = apply(ifit[ -i , ], 2, function(x) sd(x, na.rm = T)) )
    
    wdf <- as.data.frame( t(z) )
    
    actfit_long <- wdf %>% 
      mutate(
        dist = as.numeric(dist),
        drop = as.numeric(drop),
        nochange = as.numeric(nochange),
        add = as.numeric(add) ,
        utility_ego=as.numeric(utility_ego),
        utility_alter_mean = as.numeric(utility_alter_mean),
        utility_alter_sd = as.numeric(utility_alter_sd)
      ) %>%
      pivot_longer(cols = c( drop:add, utility_ego:utility_alter_sd )) 
     
    
    plt <- actfit_long %>% ggplot(aes(x=value))+ 
      # geom_density() + 
      geom_histogram() +
      facet_grid(dist ~ name, scales='free_x') + theme_bw() + 
      ggtitle(sprintf('Utility Transition Paths: Actor %s, Step %s', i, step_id))
    
    pltlist[[ length(pltlist)+1 ]] <- plt
    
    # z <- apply(iland, 1, function(x) {
    #   x <- plyr::count( x - bi_env_mat_step[i, ]) 
    #   cnts <- merge(x$freq, data.frame(x=c())
    #   names(cnts) <- x$x
    #   return(cnts)
    # }, simplify = T)
    # 
    # stringdist::stringdist( , , method='hamming')
    
    step_actor_key <- sprintf('%s|%s',step_id,i)
    
    statsl[[step_actor_key]] <- list(
      actfit_long = actfit_long,
      act_counterfacts = act_counterfacts,
      plt = plt,
      actor_id=i, chain_step_id = step_id,
      strategy=as.factor(self$strat_1_coCovar[i])
    )
    
    utilist[[step_actor_key]] <- actfit_long %>% 
      filter(name %in% c('utility_ego','utility_alter_mean')) %>%
      mutate(actor_id=i, chain_step_id = step_id, strategy=as.factor(self$strat_1_coCovar[i]))
    
  } ##/end actor loop
  
  # # Arrange plots across multiple pages
  # multi_page_plot <- marrangeGrob(pltlist, nrow = 1, ncol = 1)
  # 
  # # Save to a PDF
  # ggsave("multi_page_plots.pdf", multi_page_plot, width = 8, height = 8)
  
} ##/end step loop
## 
## step 1: actors  1  2  3  4  5  6 
## step 7: actors  1  2  3  4  5  6 
## step 13: actors  1  2  3  4  5  6 
## step 19: actors  1  2  3  4  5  6 
## step 25: actors  1  2  3  4  5  6 
## step 31: actors  1  2  3  4  5  6
#
utildf <- data.table::rbindlist(utilist, use.names = T, idcol = 'step_actor_key')


##--------------------------------------------------------------
## 4.  Actor Utility Transition Paths
##--------------------------------------------------------------
# Arrange plots across multiple pages
plt_act_paths <- marrangeGrob(pltlist, nrow = 1, ncol = 1)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
plt_act_paths

# # Save to a PDF
# actor_path_plotfile <- sprintf('%s_actor_pd_dist2peaks.pdf', self$UUID)
# plot_dir <- getwd()
# ggsave(filename = file.path(plot_dir, actor_path_plotfile), plt_act_paths, 
#        width = 8, height = 8)
env3$search_rsiena_plot_actor_utility_strategy_summary()